Constructing classical field for a Bose-Einstein condensate in arbitrary trapping 
potential; quadrupole oscillations at nonzero temperatures 
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We optimize the classical field approximation of tlie version described in J. Pliys. B 40, Rl 
(2007) for tfie oscillations of a Bose gas trapped in a harmonic potential at nonzero temperatures, 
as experimentally investigated by Jin et al. [Phys. Rev. Lett. 78, 764 (1997)]. Similarly to exper- 
iment, the system response to external perturbations strongly depends on the initial temperature 
and on the symmetry of perturbation. While for lower temperatures the thermal cloud follows 
the condensed part, for higher temperatures the thermal atoms oscillate rather with their natural 
frequency, whereas the condensate exhibits a frequency shift toward the thermal cloud frequency 
(m = mode), or in the opposite direction (m = 2 mode). In the latter case, for temperatures ap- 
proaching critical, we find that the condensate begins to oscillate with the frequency of the thermal 
atoms, as in the m = mode. A broad range of frequencies of the perturbing potential is considered. 



I. INTRODUCTION 

Experiments with atomic Bose-Einstein condensates 
driven by an external perturbation have allowed veri- 
fication of the mean-field description of the condensed 
phase in a dynamical regime where the system responds 
by collective motion. Typically, a periodic perturbation 
(with a particular symmetry) of a trapping potential was 
used to excite the Bose gas P, 2, 3, 4,], although other 
kinds of trap distortion leading, for instance, to scissors 
mode excitations 0] or the transverse monopole mode 
excitation in an elongated condensate^ were also tried. 
The investigation of the low energy collective modes of 
the condensate in the zero-temperature limit [l|, 0] has 
revealed that the mean-field description of that system 
(i.e., based on the Gross-Pitaevskii equation) works well 
[7[. However, when the study of low- lying excitations 
was extended to include the measurement of frequencies 
and damping rates function of temperature [1, 0], 
it became clear that a new theory of an interacting Bose 
gas at nonzero temperatures is required. 

The JILA experiment of Ref. Q showed two effects. 
First, it was found that two collective modes with dif- 
ferent symmetries (quadrupole modes with angular mo- 
menta equal to TO = and to = 2) behave in qualitatively 
different ways. When the temperature increases, they 
exhibit a frequency shift in opposite directions. More- 
over, for the m = mode a rather sudden upward shift 
is observed, suggesting the existence of a characteristic 
temperature which is approximately 0.65 of the critical 
temperature for the corresponding ideal gas. Secondly, 
the damping of the collective oscillations turned out 
to be dramatically dependent on temperature, showing 
that the condensate modes are damped even faster than 
the noncondensed fraction while approaching the critical 
temperature. All these puzzling findings triggered a lot 



of theoretical work and after a few years resulted in the 
development of Zaremba-Nikuni-Griffin formalism [1, Q 
and the second-order quantum field theory [lol . for a 
Bose gas. Recently, another attempt to describe finite- 
temperature properties of low-lying collective modes was 
undertaken in Ref. [T^ within the approximation called 
the projected Gross-Pitaevskii equation. 

The Zaremba-Nikuni-Grifhn formalism [9| applied to 
the results of the JILA experiment gives relatively good 
agreement. Also the calculations based on the second- 
order quantum field theory [l0| show a good agreement 
with experimental data. However, already these two ap- 
proaches differ when considering their fundamentals. For 
example, the first one neglects the phonon character of 
the low-lying energy modes, the anomalous average, as 
well as the Beliaev processes. Regarding the JILA ex- 
periment, for the TO = mode, the Zaremba-Nikuni- 
Grifhn method predicts an additional branch of conden- 
sate frequencies, so far not observed experimentally. No 
such branch is found within the second-order theory of 
Ref. [l^. Furthermore, this approach predicts a single 
frequency of condensate response for a particular tem- 
perature, regardless of driving frequency, and hence the 
notion of a natural condensate frequency is meaningful. 
This strongly differs from what is reported in Jackson 
and Zaremba Q , where the condensate response depends 
on the driving frequency of the whole system (condensed 
and noncondensed parts). A comprehensive discussion of 
both approaches in the context of JILA experiment can 
be found in a recent review article p^ . Another formal- 
ism, based on the projected Gross-Pitaevskii equation, 
used to model the JILA experiment 3 produces good 
agreement with experimental data up to 0.65Tc, and for 
the TO = 2 mode at higher temperatures. However, it 
fails to predict the sudden upward frequency shift for the 
TO = mode at this temperature |12l ]. This failure is 
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perhaps related to the way the cutoff parameter (which 
splits the space of modes into the highly occupied ones 
that are described in terms of the classical field, and the 
others that are sparsely occupied and require, in princi- 
ple, a quantum treatment) is chosen. The details of the 
splitting procedure within the projected Gross-Pitaevskii 
equation approach are discussed in Ref . • 

In this paper we apply the classical field approxima- 
tion in the version described in Ref [l^ to the case of 
a trapped interacting Bose gas at nonzero temperatures 
driven by an external perturbation as in the experiment 
of Ref. Q. The main purpose of this work is to check 
whether the classical field approximation is able to re- 
produce, at least qualitatively, the findings of JILA ex- 
periment. This is an especially important task because 
of the recently reported failure |l2| to explain the be- 
havior of the m — mode within the projected Gross- 
Pitaevskii equation method, which is conceptually very 
close to our classical field approximation. On the other 
hand, although the other existing theories and [3 
both lead to relatively good agreement with the experi- 
ment , they somehow contradict each other conceptu- 
ally as discussed in the previous paragraph. It would be 
nice to have an alternative view of the processes going 
on in the perturbed Bose gas. Finally, the approaches [9| 
and [lol | have some difficulties in describing the dynamics 
of the thermal cloud, especially in the m = 2 mode. In 
fact, no predictions for thermal component frequencies 
in the m = 2 mode are given in 0] or [lol |. Simulta- 
neously, within the projected Gross-Pitaevskii equation 
method, at higher temperatures, the thermal cloud (in 
fact, for both the m = and m = 2 modes) oscillates at 
frequencies much lower than in experiment. 

The classical field approximation has already been ap- 
plied in static and dynamical regimes for a uniform and 
harmonically trapped systems (for a review, see [3). 
This approach was used to inv estig ate the thermody- 
namics of an interacting gas [igI . [l7l | as well as dynami- 
cal processes like the photoassociation of molecules [l^ , 
the dissipative dynamics of a vortex [l^, the superflu- 
idity in ring-shaped traps [20l |. or the thermalization in 
spinor condensates [2l|. The classical field approxima- 
tion was also tested at a quantitative level when e.g. 
the Bogoliubov-Popov quasiparticle energy spectrum in 
a uniform Bose gas was obtained or when the pro- 
cess of splitting of doubly quantized vortices in dilute 
Bose-Einstein condensates [22| was studied. 

The paper is organized as follows: In Sec. |TT] we de- 
scribe the classical field approximation for a trapped 
Bose gas with particular attention on how the equilib- 
rium state is obtained and on the quality of the solu- 
tion in comparison with the equilibrium states calculated 
within the self-consistent Hartrce-Fock method. Sec. IIIII 
discusses the results (response frequencies and damping 
rates) for the quadrupole m = collective mode for pa- 
rameters as in JILA experiment Q while in Sec. IIVI we 
do the same for m ~ 2 excitation. Finally, we conclude 
in Sec. El 



II. CLASSICAL FIELD APPROXIMATION FOR 
A TRAPPED BOSE GAS 

A. Formalism 

A good starting point to introduce the classical field 
approximation is the usual Heisenberg equation of mo- 
tion for the bosonic field operator *I'(r,i) which annihi- 
lates an atom at point r and time t. The field operator 
^(r,t) fulfills standard commutation relations: 



*(r,<),*+(r',i) =(5(r-r') 



(1) 



with other equal time commutation relations for [5',^'] 
and [^'+, ^P+J being zeros. The equation of motion reads: 



d ~ 
ih—-^{Y,t) 



-^V^ + F.(r,t) 



*(r,i) 



-.9*+(r,t)^'(r, t)^{Y, t). 



(2) 



Here, we assume the time-dependent trapping poten- 
tial Vtr and the usual contact interaction for colliding 
atoms. The coupling constant g = ATrh^a/m is expressed 
in terms of the s-wave scattering length a. 

Next, we expand the field operator ^(r,t) in the basis 
of one-particle wave functions '(/'^(r), where fc is a set of 
one-particle quantum numbers: 



(3) 



Now, we assume that some of the modes used in the ex- 
pansion ^ are macroscopically occupied and extend the 
original Bogoliubov idea [1^ by replacing all operators 
ak{t) corresponding to these modes by c-numbers. When 
only macroscopically occupied modes are considered, the 
field operator ^(r,t) is turned into the complex wave 
hmction 5'(r,t) and the expansion ([3]) takes the form: 



(4) 



The upper index in the summation tells us that, indeed, 
the wave function ^'(r, t) is expanded only over a finite 
number of states, i.e. those which are macroscopically 
occupied. We call the wave function 5'(r, i) the classical 
field. This is analogous to the way that intense electro- 
magnetic waves can be treated. In spite of consisting of 
photons, an intense light beam is well characterized by 
the classical electric and magnetic fields. Since the ex- 
periments with dilute atomic gases are performed with 
millions of atoms it seems to be a plausible approxima- 
tion to use the classical field to describe atoms in analogy 
with electric and magnetic fields for photons. 

Obviously, the classical field obeys the following equa- 
tion: 



.4v,(r,t) = 



-|^V^ + F.(r,t) 



*(r,t) 



+.g**(r, t)*(r,t)^'(r, t) 



(5) 
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We usually implement the cutoff parameter kmax by solv- 
ing the Eq. ([5]) on a rectangular grid using the Fast 
Fourier Transform technique. The spatial grid step de- 
termines the maximal momentum per particle (and hence 
the energy) in the system whereas the use of the Fourier 
transform implies a projection in momentum space. 

The equation (O looks like the usual Gross-Pitacvskii 
equation describing the Bosc-Einstcin condensate at zero 
temperature. However, here the interpretation of the 
complex wave function 5'(r,i) is different. It describes 
all the atoms in the system, both those in a condensate 
and in a thermal cloud. Therefore, the question appears 
how to split the classical field ^'(r, t) into the condensed 
and noncondensed fractions. For that we use the defini- 
tion of a Bosc-Einstein condensation proposed originally 
by Penrose and Onsager [23| . According to this definition 
the condensate is assigned to the eigenvector correspond- 
ing to the dominant eigenvalue of a one-particle density 
matrix. So, we built the one-particle density matrix in 
the following way: 

p«(r,r';t) = lvl/(r,t)vl/*(r',i), (6) 

where N is the total number of particles. However, and 
here comes the surprise ~- since ([6]) is just the spectral de- 
composition of a one-particle density matrix, this would 
imply that the classical field ^'(r, i) is the condensate 
wave function and that all atoms are in the condensate. 
To split the classical field into the condensed and thermal 
fractions one first realizes that the high energy solutions 
of Eq. ^ oscillate rapidly in time and space. On the 
other hand, the detection process is always performed 
with limited spatial and temporal resolution. Therefore, 
it becomes clear that the measurement process with its 
limited resolution involves a kind of averaging (coarse- 
graining) of Eq. ([6]). Again, an analogy with electro- 
magnetic waves is in place. A partially coherent light 
beam exhibits complicated spatial and temporal behav- 
ior on short scales, in fact too complicated to be typically 
measured. What is important are then the correlation 
functions at long enough spatial and temporal separa- 
tion, averaged over the smaller time intervals and space 
segments. Similarly here, the averaging over time and/or 
space of the one-particle density matrix ([6]) results in a 
partial loss of the information contained in the classical 
field [TB|. In other words, the mixed state emerges out of 
the pure one. 

In a typical experiment, what is measured is the col- 
umn density along some direction. Hence, we also imple- 
ment the coarse-graining procedure in our numerics in 
this manner as: 

p{x, y, x', y'; J *(a;, y, z, t) **(a;', y', z, t) . 

. 

Solving the eigenvalue problem for a coarse-grained 
density matrix ([7]) results in a decomposition: 

P = ^nkipk{x,y,t)ipl{x\y',t) , (8) 

k 



where Uk = Nk /N are the relative occupations of macro- 
scopically occupied modes i-pk- Defining the dominant 
eigenvalue as no, the condensate wave function (normal- 
ized to No/N) can be written as: 

r/Vo 

'>Po{x,y,t) = ]l — (po{x,y,t) . (9) 

All the other modes contribute to the thermal density 
which is, therefore, given by: 

PT{x,y,t) = p{x,y,x,y;t) - \ipo{x,y,t)f . (10) 

The appropriateness of the averaging ^ for obtain- 
ing the large eigenvalues of a coarse-grained density ma- 
trix was verified by us in various ways. For example, 
we checked that for the classical field at equilibrium this 
procedure gives the same results as averaging over a long 
enough time. Prescription ([7]) also works at zero temper- 
ature (when all atoms are expected to be in a conden- 
sate), having been successfully tested in this respect in 
the demanding case of a lattice of bent vortices (see Ref. 
[25! for further details). 

According to the description given above the splitting 
of the system into condensed and noncondensed compo- 
nents is a result of Bose statistics, interaction, and the 
measurement process. Unlike the alternative approaches 
@j we do not impose a two-component character of 
the system from the beginning. Also, our version of the 
classical field method is well suited to describe single re- 
alizations of the experiment since it corresponds to a 
microcanonical ensemble, as opposed to the competing 
approach (26| . which deals with canonical ensembles. 

B. Obtaining equilibrium states 

An initial classical field is generated from the ground 
state solution of Eq. ([5|) by adding appropriate random 
disturbance. An inspection of that equation indicates 
that only the product of the coupling constant g and the 
total number of atoms N enters it. Therefore, we nor- 
malize the initial classical field 5'(r, i = 0) to unity. The 
norm of \E'(r, t) is then one of the constants of motion of 
Eq. ((S]). Another constant of motion is the total energy. 
Such an initial state is then propagated according to Eq. 
([5]) imtil the constituent energies (kinetic, trap, and inter- 
action) cease to change systematically in time and exhibit 
only fluctuations. In this way the classical field at ther- 
mal equilibrium, corresponding to the particular values of 
gN and Etot/N, is obtained. An example is given in Fig. 
[U where we plot cuts of the total density, the condensate, 
and the thermal densities according to the prescription 
detailed in the previous section. Here, the equilibrium 
classical field describes a degenerate ^^Rb Bose gas in a 
magnetic trap with frequencies lo± = i^x,y = 27r x 129 Hz 
and Lu^ = 2tt X 365 Hz as in the experiment of Ref. 0] . 
Other parameters are: gN = 2911.9 and Etot/N = 21.2 
in units of hujzlh/mujz)^^'^ and hwz, respectively. They 
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result in a condensate fraction no = 0.3. A bimodal dis- 
tribution (thick solid line) is clearly visible. 
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FIG. 1: Cuts of the total (thick solid line), condensate (dashed 
line) , and thermal (thin solid line) densities as obtained by in- 
tegrating a one-particle density matrix along the direction of 
imaging and then applying the splitting procedure described 
in Sec. Ill Al Note the bimodal character of the density dis- 
tribution (here, the condensate fraction equals 0.3). The os- 
cillatory unit of length is defined via the axial trap frequency: 
— and equals 0.565 fj,m. 



Now we have to solve the problem how to find the 
number of particles assigned to the classical field at equi- 
librium and how to determine the temperature of the 
system. This can be done in two ways. The first is 
given here, the second in Sec. Ill Cl To begin with, having 
the classical field at equilibrium one can project the field 
^'(r, t) on the harmonic oscillator states obtaining in this 
way the relative populations of these states. Fig. [2]shows 
relative populations for gN = 1811.9 and Etot/N = 10.0 
as a function of the harmonic oscillator states' energy 
(precisely, the time average over 274 ms is plotted). The 
maximal one-particle mode energy is determined by the 
momentum cutoff Pmax as ^^^^./m. An important ob- 
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FIG. 2: Relative populations of various harmonic oscillator 
states. The energy cutoff equals Pmax/m, where Pmax is the 
momentum cutoff. 



servation is made when one looks at the energy accu- 
mulated in harmonic oscillator states, see Fig. [3l For 
higher energy states the product Si (where and Si 
are the relative population and the harmonic oscillator 
state energy, respectively) becomes constant. On the 
other hand, for highly occupied modes (i.e., modes sat- 
isfying £i — ^J^ ksT, where fj. is the chemical potential) 
the quantum Bose-Einstein distribution reduces to the 
classical equipartition distribution. For the classical field 
studied here the equipartition extends all the way to the 
cutoff energy: 



{ei - Ai) = ksT , 
which can be written equivalently as 

Hi (e,; - fJ.) = ksT/N . 



(11) 



(12) 



Therefore, since energy equipartition is established, the 
higher energy harmonic oscillator states become the 
quasiparticle modes. In Fig. [3] we determine the ratio 
ksT/N to be 4.62 x 10"''. Having the ratio ksT/N we 
now find from Eq. \12\ the relative populations of high 
energy quasiparticle modes. In particular, we get the 
relative population of the least occupied modes which 
belong to the classical field. In the example here, they 
are 2.99 x 10"^ (based on formula ^ with ksT/N ob- 
tained above) . In Ref . [l^l , arguments are given that the 
best match between this method and the ideal Bose gas 
occurs when the occupation of the least occupied mode is 
Ncut = 0.46. These are based on the comparison between 
the probability distribution of the ideal Bose gas and its 
classical field counterpart. Assuming, then that the aver- 
age number of atoms in this least-occupied mode is 0.46 
we can retrieve the total number of atoms in the system 
and its temperature separately. For the example here, 
these are TV = 15342 and T = 124.1 nK. Since no data 
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FIG. 3: Relative populations multiplied by the state ener- 
gies for various harmonic oscillator states. The figure clearly 
shows that for higher energy states the equipartition of an 
energy is established. Therefore, the higher energy states be- 
come quasiparticle modes. The fraction ksT/N equals ap- 
proximately 4.62 X 10~*. 



on the population of the cutoff mode are available for a 
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weakly interacting Bose gas considered here we compare 
in Table H] the temperature of the system and the total 
number of atoms for various values of Ncut- 



N^ut 


T[nK] 


TV 


0.44 


119.5 


14778 


0.46 


124.1 


15342 


0.48 


130.4 


16121 


0.50 


135.8 


16793 


0.52 


141.2 


17465 



TABLE I: Temperature and the total number of atoms ob- 
tained by projecting the classical field on the harmonic oscil- 
lator states of the example described in Sec. Ill Bl for different 
values of cutoff parameter. 



Comparison with the self-consistent 
Hartree-Fock model 



Another approach to obtain the number of atoms and 
the temperature of the system described by the classi- 
cal field at equilibrium is based on the self-consistent 
Hartree-Fock method [1^ . Since this method works well 
for a Bose gas at equilibrium as verified experimentally in 
[2^, a comparison between this model and the classical 
field approximation should be instructive. The Hartree- 
Fock description is defined by the following set of equa- 
tions nil: 

nc{r) = -[^i- Vtr{r) - 2.gnt^(r)] (13) 
9 

/(r, p) = (e[pV2™+v.„(r)-M]/fei.T _ -1 (^4) 



nth{r) = 33/2 (e 



[^L~v,ff{v)]/kBT 



Veffir) = Vtrir)+2gncir) + 2gnthir) 
^ = gndO) + 2gnthi0) , 



where 



At = 



rnksT 



1/2 



(15) 

(16) 
(17) 



(18) 



is the thermal de Broglie wavelength and the 33/2 (^) 
function is given by the expansion: 



(Us])). Since both the effective potential, 14//(r), and the 
chemical potential, fi, appearing on the right-hand sides 
of Eqs. (fO)l . and depend on the condensate 

and thermal densities the set of Eqs. and is 

well suited to be solved iteratively. For that, however, 
we have to first choose the temperature of the system 
and then keep fixed the total number of atoms (which is 
N — J dr {nc{r) + nth{y^))) during the iterations. Hav- 
ing the densities of condensed and noncondensed frac- 
tions one can easily calculate two important parameters: 
the condensate fraction and the total energy per particle. 
Now the strategy is as follows: find the input parameters 
{N, T) in such a way that the condensate fraction and the 
total energy per atom calculated within the Hartree-Fock 
method match the values calculated from the classical 
field at equilibrium. This procedure allows one to deter- 
mine the number of atoms and the temperature assigned 
to the classical field separately. 

These parameters for the example discussed in the con- 
text of Figs. [2] and [3] are found to be iV = 17306 and 
T = 128.7 nK. These values are very close to what was 
obtained in the previous section with the cutoff occupa- 
tion Ncut = 0.46, and differ by 3.7% in the temperature 
and 12.8% in the total number of atoms. The agreement 
is good even though the average occupation of the highest 
energy modes (the last modes considered in the classical 
field approximation as being macroscopically occupied) 
was taken the same as for an ideal gas. Unfortunately, 
there is no data available for the average occupation of 
the cutoff region modes for the weakly interacting gas 
considered here. Changing slightly this cutoff occupa- 
tion number, the agreement between both approaches to 
the system parameters can be made even better (for ex- 
ample, for Ncut = 0.48 the difference is 1.3% and 7.4% in 
the temperature and the total number of atoms, respec- 
tively). In Fig. m the total Hartree-Fock and classical 
field densities are plotted together for parameters as in 
Figs. [2] and [3] showing a good agreement. The classical 
field density, however, exhibits all the fluctuations which 
arc not present in the Hartree-Fock model. 

In Table |ll] we again compare the parameters of the 
system (total number of atoms and temperature) for the 
methods described in this and in the previous sections. 
However, this time the comparison is for several different 
equilibrium states. All the states being compared are 
used later in the simulations related to JILA experiment 
(3 (see Sees. IH] and |IV| . 



ff3/2(z) X! 



,3/2 



(19) 



The basic variables in this approach are the conden- 
sate density, nc(r), and the distribution function in phase 
space for thermal atoms, /(r,p). They are calculated 
according to the Eqs. P^ and (HH). The thermal den- 
sity, ntft(r), is just an integral of the distribution function 
/(r, p) over momenta and can be found analytically (Eq. 



D. Obtaining equilibrium states on demand 



In Sees. Ill Bl and III CI we detailed the methods for re- 
trieving the total number of atoms and the temperature 
assigned to the classical field at equilibrium. Since the 
product gN is an initial parameter it means that the cou- 
pling constant g (and consequently the scattering length 
a) is known only after the classical field is thermalized. 
Unfortunately, it usually differs from the value that was 
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X (osc. units) 

FIG. 4: Total density cuts for the self-consistent Hartree- 
Fock (thick line) and classical field (thin line) methods for a 
system with the total number of atoms = 17306 and the 
temperature T = 128.7 nK. 





CPA 


HF 


r[nK] 
N 


79.8 
6751 


91.8 
11300 


T[nK] 
N 


124.1 

15342 


128.7 
17306 


T[nK] 
N 


142.1 
16347 


153.9 
18699 


T[nK] 
N 


154.4 
19313 


168.8 
21509 


T[nK] 
N 


177.8 

26984 


190.4 
27875 


T[nK] 
N 


195.3 
27770 


204.8 
31222 


T[nK] 
N 


264.1 
61876 


240.9 
43572 



TABLE II: Comparison (for several final equilibrium states as 
per Table lln|l between the results (the temperature and the 
total number of atoms) obtained by projecting the classical 
field on the harmonic oscillator states as described in Sec. IIIBI 
(CFA column) and by utilizing the self-consistent Hartree- 
Fock method according to Sec. Ill CI (HE column). For CFA 
data Ncut = 0.46. 



used to calculate an initial value of gN . Therefore, an ap- 
proach for obtaining the classical field at equilibrium cor- 
responding to given values of the total number of atoms 
and the temperature is required. 

The method we have developed is rather simple al- 
though demanding from a numerical point of view. Let's 
assume that we need the classical field which at equi- 
librium describes the system with given parameters TV 
and T. We start from a solution obtained from the 
self-consistent Hartrce-Fock method corresponding to the 
chosen parameters. Then we build an initial classical field 



as follows: 

^{v,t = Q) = ^J^) + ./^^)e'''^'''^ (20) 

and randomize the phase ^p{v) and the density rithir) in 
such a way that the total energy per atom in the classical 
field equals the corresponding energy in the Hartree-Fock 
model. The presence of the phase factor in the second 
term in pO|) is necessary. Without this, the classical field 
suffers from a lack of kinetic energy in comparison with 
the Hartree-Fock method, where it is calculated from the 
distribution fimction /(r,p): 

Ekrn = -^ J drdp^f{r,p) 

= ^-^ j '''''' A''"^'"'^'"")- 

Now we evolve the classical field according to Eq. ([5|) 
and let the field thermalize. During the thermalisation, 
the total energy per atom is a constant of motion, but 
the condensate fraction usually changes. However, there 
will be a particular value of the spatial step of the grid for 
which the condensate fraction does not change in time. 
Then, since the total energy per atom is a constant of mo- 
tion the values of parameters tiq and Etot /N at the end of 
thcrmalization process are the same as at the beginning. 
Consequently, the number of particles N and the tem- 
perature T must be the same as chosen at the beginning 
as well. Although, the procedure just described is nu- 
merically time consuming (since it requires several trials 
to obtain the proper value of the spatial step), it is much 
more efficient than attempts to match final T and N si- 
multaneously, which require search in two-dimensional 
parameter space. Here, the energy matching is computa- 
tionally fast because it requires no thermalisation, while 
the final T, N matching is done with only one free pa- 
rameter, the spatial lattice spacing. 

III. RESULTS FOR THE M=0 MODE 

Our numerical procedure involves the following steps: 
First, we find the classical field (as described in lH Pp cor- 
responding to the Bose gas at equilibrium confined in a 
harmonic trap with frequencies = i^x,y = 27r x 129 Hz 
and ujz = 2t: X 365 Hz. According to the experiment Q 
the total number of atoms was of the order of ten to a few 
tens of thousands and the initial temperature was ranging 
up to the critical temperature. The number of condensed 
atoms remained on a level of several thousand. The nu- 
merical values of A^, Nq, T, and the reduced temperature 
T' = T/Tc (with Tc being the transition temperature for 
a harmonically confined ideal gas) for the states used in 
the simulations are shown in Table [Till 

Next, the Bose gas is disturbed by a sinusoidal pertur- 
bation to the trapping potential. Since the classical field 
describes both the condensed and noncondensed atoms, 
the disturbance of the classical field means that both 
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N 


No 


T[nKl 


T' 


11300 


8558 


92.8 


0.497 


17306 


10349 


128.7 


0.605 


18700 


7905 


153.9 


0.705 


21509 


7768 


168.8 


0.738 


27875 


8496 


190.4 


0.763 


31222 


7753 


204.8 


0.791 


43572 


7111 


240.9 


0.832 



TABLE III: Numerical values of iV, No, T, and T' used in 
the simulations. 

components are simultaneously perturbed. To excite the 
TO = and the to = 2 quadrupole modes, the perturba- 
tion of the trapping potential takes the form: 

SVtr (r , t) A{t) [luIx^ sin(wd t + (/))+ Lo^y^ sin(wd t)] , 

(22) 

where Ud is the driving frequency and </> is a phase shift 
between the x and y direction perturbations. The choice 
of the phase shift (j) determines the symmetry of the ex- 
cited collective mode - for ^ = (0 = tt) the to = 
{m = 2) mode is excited. The perturbation, as in the ex- 
periment, lasts for 14 ms and the amplitude A{t) (= 0.05) 
takes small value to avoid any nonlinear effects. 

After the perturbation is turned off, the classical field 
is oscillating in time. Is it possible that the condensed 
and noncondensed components (extracted from the single 
classical field) exhibit oscillations with different frequen- 
cies? To answer these question we split the classical field 
into the condensate and the thermal cloud in the way 
described in Sec. Ill Al and calculate the widths of both 
components from the formula: 

Wc,th = J dxdy {x^ + y'^)nc,th{x,y) ■ (23) 

Results are shown in Fig. [S]for the to = mode. Solid 
symbols (upper frame) represent the condensate widths, 
whereas the open symbols (middle frame) stand for ther- 
mal cloud widths. Solid lines in both frames are fits by 
an exponentially damped sine waves: 

Ac^th exp(-7c,tht) sin(wc,t/ii -I- -dc^th) + Beth (24) 

to numerical data. As in the experiment, fits are per- 
formed based on five initial oscillations. The lower frame 
shows, indeed, that condensate and thermal components 
oscillate with different frequencies. Moreover, these os- 
cillations are phase shifted and the condensate oscillates 
slower than the thermal cloud. 

In Figs. [6] and [7] wc summarize our results for the 
TO = mode. In Fig. [B]we plot the frequencies of the 
condensate and thermal fractions response to the exter- 
nal perturbation for various temperatures together with 
the experimental data of Ref. Q ■ Black solid and open 
symbols represent data for a condensate (upper frame) 




TIME (ms) 

FIG. 5: Widths of a condensate (upper frame) and its ther- 
mal cloud (middle frame) as a function of time for the m = 
mode. The widths (solid and open symbols) are calcu- 
lated from the formula (|23p after the condensate and the 
thermal cloud are extracted out of the classical field. The 
solid lines are fits to an exponentially damped sine function: 
Ac^th expi-Jctht) sin{ijjc,tht + ^c,th) + Bc,th- These fits al- 
low us to obtain the frequencies and the damping rates of 
the oscillations of the condensate and the thermal cloud in 
response to the external perturbation. The initial temper- 
ature of the cloud and the driving frequency are T' — 0.8 
and ujd ~ 1 .75 Lu±. The lower frame shows more clearly the 
frequency and phase shifts between the condensate and the 
thermal cloud. 



and a thermal cloud (lower frame) for various driving 
frequencies according to the legend attached to the fig- 
ure. Gray symbols with error bars are the experimental 
results. Up to temperature T w O.BTc both components 
oscillate with the same frequency, which is the natural 
condensate frequency for the to = collective mode. At 
approximately 0.65Tc a rather sudden upward shift in 
condensate frequency is observed in the experiment. Our 
numerics shows that at this temperature the dynamics 
of the thermal cloud changes. The thermal component 
starts to oscillate with a higher frequency approaching 
eventually 2uj± which is the oscillation frequency of a 
thermal gas alone. For higher temperatures the thermal 
fraction becomes dominant and finally thermal atoms 
change the dynamics of the condensate in such a way that 
condensed atoms oscillate along with the thermal ones. 
Unfortunately, the influence of the thermal cloud on the 
condensate is apparently too weak and consequently the 
condensate starts to oscillate with higher frequencies only 
for higher temperatures. 

Our results also show that the notion of natural con- 
densate frequency breaks down when the thermal cloud 
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FIG. 6: Condensate (upper frame) and the thermal cloud 
(lower frame) frequencies for the m = collective mode as a 
function of reduced temperature. Black (solid and open) sym- 
bols correspond to the numerical results obtained for various 
driving frequencies (according to the legend) whereas gray 
symbols (solid for the condensate and open for the thermal 
cloud) with error bars are taken from the experiment by Jin 
et al. 01 . 



FIG. 7; Damping rates for the condensate (upper frame) and 
the thermal cloud (lower frame) for the m = mode as a func- 
tion of reduced temperature. Black (solid and open) symbols 
correspond to the numerical results obtained for various driv- 
ing frequencies (according to the legend) whereas gray sym- 
bols (solid for the condensate and open for the thermal cloud) 
with error bars are taken from the experiment by Jin et al. 
a- 



is present. The condensate response depends on the dy- 
namics of the thermal cloud and. in fact, the possible fre- 
quencies for the condensate oscillation lie in an interval 
which gets wider for higher temperatures. In particular, 
no two branches of frequencies are visible as reported in 
Ref. i. In Fig. [7] we compare the numerical and ex- 
perimental damping rates for the oscillations of the con- 
densed and noncondensed components. There is some 
discrepancy for higher temperatures where the thermal 
component seems to be damped too strongly in compar- 
ison with the experimental data. Perhaps this increases 
the temperature at which the frequency of the condensate 
fraction shifts up. 



IV. RESULTS FOR THE M=2 MODE 

We now do the same analysis for the m = 2 mode. The 
system is disturbed by changing the trapping potential 
according to (|22p with the phase shift <j> — it. Such a 
perturbation excites the quadrupole m = 2 oscillations. 
The classical field is split into the condensed and ther- 
mal parts and the condensate and thermal widths are 
calculated using formulas exhibiting the m ~ 2 mode's 
symmetry: 

(25, 

As before, these data are fitted by exponentially damped 
sine waves. 

Frequencies and damping rates for the condensate and 
the thermal cloud are plotted in Figs. [5] and [HI Again, 



the system responds in a way that depends on the driving 
frequency. A comment on how the authors of the exper- 
iment [3| choose the driving frequency is in place here. 
They say that the driving frequency is set to match the 
frequency of the excitation being studied. This could 
make sense for m = mode where supposedly the con- 
densed and noncondensed fractions always oscillate in 
phase. Supposedly, because in Q we see only two data 
points (both having frequency approximately 2uj±) show- 
ing the behavior of the thermal cloud in the m — 
mode. However, for the m = 2 mode the frequencies 
at which the condensate and thermal cloud oscillate are 
different, and the meaning of the driving frequency as 
the frequency of the excitation being studied is not clear. 
Therefore, in Fig. [8] we show our data for various driving 
frequencies. Good agreement with experiment is found 
for appropriate driving frequencies up to temperatures 
~ O.STc. For temperatures approaching the critical one, 
however, we observe the same effect as for the m = 
mode. The thermal cloud becomes dominant and the 
frequencies at which the condensate oscillates become 
higher. Fig. [8] shows that already for a temperature 
T = 0.83rc, with the driving frequency ujd = l.75uj±, 
the condensate oscillates with the frequency of the ther- 
mal cloud. Such a behavior was not observed experimen- 
tally. On the other hand, it is an expected behavior, since 
the condensed fraction becomes smaller and smaller while 
the critical temperature is approached, and the dynamics 
should be dominated by the thermal cloud. In the case 
of the thermal cloud, our model predicts frequencies in 
agreement with the experiment (other theoretical studies 
of the JILA experiment either do not show results for the 
thermal atoms for the m = 2 mode, or predict frequencies 
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that are very different than observed experimentally). 

In Fig. [9]wc compare the numerical and experimental 
damping rates for the oscillations of the condensed and 
noncondensed parts of the system. As for the m = 
mode, there is some discrepancy for higher temperatures 
where the thermal cloud is damped too strongly, whereas 
the condensate is damped too weakly. 
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FIG. 8: Condensate (solid symbols) and thermal cloud (open 
symbols) frequencies for the m = 2 collective mode, as a func- 
tion of reduced temperature. Black symbols correspond to the 
numerical results obtained for various driving frequencies (ac- 
cording to the legend, here the shape of the symbol determines 
the driving frequency), whereas gray symbols with error bars 
are taken from the experiment by Jin et al. [sj. Note, how- 
ever, that the only experimental point for the thermal cloud 
is marked with an open gray circle. 



CONCLUSIONS 



In conclusion, we have presented in detail the construc- 
tion of the classical field describing the desired number 
of atoms, confined in any trapping potential, at a pre- 
scribed temperature. We have studied the oscillations of 
the Bose-Einstein condensate in the presence of a ther- 
mal cloud. As in the experiment, we find the tempera- 
ture dependent condensate frequency shift for both the 
m = and m = 2 collective oscillation modes. For the 
m = mode, the thermal atoms pull the condensate 
fraction along, and above some characteristic tempera- 
ture (« Q.STc) the condensate tends to oscillate with the 
frequency of the thermal part (approximately equal to 
2w_l). Unfortunately, in the present version of the clas- 
sical field approximation, the value of this characteristic 
temperature turns out to be about 20% higher than the 
one observed in the experiment. For the m = 2 mode, on 
the other hand, the frequency at which the condensate 
oscillates first decreases with an increase of temperature 
(the thermal cloud oscillating with its natural frequency 
equal to 2uj±_ damps the condensate motion) but when 
the temperature gets closer to the critical temperature 
the condensate starts to oscillate with higher frequencies, 
approaching the frequency of the thermal cloud. 
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FIG. 9: Damping rates for the condensate (upper frame) and 
the thermal cloud (lower frame) for them = 2 mode, as a 
function of reduced temperature. Black (solid and open) sym- 
bols correspond to the numerical results obtained for various 
driving frequencies (according to the legend), whereas gray 
symbols (solid for the condensate and open for the thermal 
cloud) with error bars are taken from the experiment by Jin et 
al. Note that the only experimental point for the thermal 
cloud is marked with an open gray circle. 
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